Covariance functions: Exponential family, radial basis functions,cubic spline, compactly supported Wendland family and stationary covariances.

Description

Given two sets of locations these functions compute the cross covariance matrix for some covariance families. In addition these functions can take advantage of spareness, implement more efficient multiplcation of the cross covariance by a vector or matrix and also return a marginal variance to be consistent with calls by the Krig function.

Note: These functions have been been renamed from the previous fields functions using 'Exp' in place of 'exp' to avoid conflict with the generic exponential function (exp(...))in R.

Usage

Exp.cov(x1, x2, theta = rep(1, ncol(x1)), p = 1, C = NA, marginal=FALSE)
Exp.simple.cov(x1, x2, theta =1, C=NA,marginal=FALSE)
Rad.cov(x1, x2, p = 1, m=NA, with.log = TRUE, with.constant = TRUE, 
               C=NA,marginal=FALSE, derivative=0)
cubic.cov(x1, x2, theta = 1, C=NA, marginal=FALSE) 
Rad.simple.cov(x1, x2, p=1, with.log = TRUE, with.constant = TRUE, 
               C = NA, marginal=FALSE)
stationary.cov(x1, x2, Covariance="Exponential", Distance="rdist",
               Dist.args=NULL, theta=1.0,V=NULL,C=NA, marginal=FALSE,derivative = 0,...)
stationary.taper.cov(x1, x2, Covariance="Exponential", 
           Taper="Wendland", 
           Dist.args=NULL, Taper.args=NULL, 
           theta=1.0,V=NULL, C=NA, marginal=FALSE,
           spam.format=TRUE,verbose=FALSE,...)
wendland.cov(x1, x2, theta = 1, V=NULL, k = 2, C = NA, 
             marginal =FALSE,Dist.args = list(method = "euclidean"), 
             spam.format = TRUE, derivative = 0, verbose=FALSE)

Arguments

Value

If the argument C is NULL the cross covariance matrix is returned. In general if nrow(x1)=m and nrow(x2)=n then the returned matrix will be mXn. Moreover, if x1 is equal to x2 then this is the covariance matrix for this set of locations.

If C is a vector of length n, then returned value is the multiplication of the cross covariance matrix with this vector.

Details

For purposes of illustration, the function Exp.cov.simple is provided in fields as a simple example and implements the R code discussed below. List this function out as a way to see the standard set of arguments that fields uses to define a covariance function. It can also serve as a template for creating new covariance functions for the Krig and mKrig functions. Also see the higher level function stationary.cov to mix and match different covariance shapes and distance functions.

A common scaling for stationary covariances: If x1 and x2 are matrices where nrow(x1)=m and nrow(x2)=n then this function will return a mXn matrix where the (i,j) element is the covariance between the locations x1[i,] and x2[j,]. The exponential covariance function is computed as exp( -(D.ij)) where D.ij is a distance between x1[i,] and x2[j,] but having first been scaled by theta. Specifically if theta is a matrix to represent a linear transformation of the coordinates, then let u= x1%*% t(solve( theta)) and v= x2%*% t(solve(theta)). Form the mXn distance matrix with elements:

D[i,j] = sqrt( sum( ( u[i,] - v[j,])**2 ) ).

and the cross covariance matrix is found by exp(-D). The tapered form (ignoring scaling parameters) is a matrix with i,j entry exp(-D[i,j])*T(D[i,j]). With T being a positive definite tapering function that is also assumed to be zero beyond 1.

Note that if theta is a scalar then this defines an isotropic covariance function and the functional form is essentially exp(-D/theta).

Implementation: The function r.dist is a useful FIELDS function that finds the cross Euclidean distance matrix (D defined above) for two sets of locations. Thus in compact R code we have

exp(-rdist(u, v))

Note that this function must also support two other kinds of calls:

If marginal is TRUE then just the diagonal elements are returned (in R code diag( exp(-rdist(u,u)) )).

If C is passed then the returned value is exp(-rdist(u, v)) %*% C.

Some details on particular covariance functions:

[object Object],[object Object],[object Object]

About the FORTRAN: The actual function Exp.cov and Rad.cov call FORTRAN to make the evaluation more efficient this is especially important when the C argument is supplied. So unfortunately the actual production code in Exp.cov is not as crisp as the R code sketched above. See Rad.simple.cov for a R coding of the radial basis functions.

Examples

Run this code

# exponential covariance matrix ( marginal variance =1) for the ozone
#locations 
out<- Exp.cov( ozone$x, theta=100)

# out is a 20X20 matrix

out2<- Exp.cov( ozone$x[6:20,],ozone$x[1:2,], theta=100)
# out2 is 15X2 matrix 

# Kriging fit where the nugget variance is found by GCV 
# Matern covariance shape with range of 100.
# 

fit<- Krig( ozone$x, ozone$y, Covariance="Matern", theta=100,smoothness=2)

data( ozone2)
x<- ozone2$lon.lat
y<- ozone2$y[16,]
# Omit the NAs
good<- !is.na( y)
x<- x[good,]
y<- y[good]


# example of calling the taper version directly 
# Note that default covariance is exponential and default taper is 
# Wendland (k=2).

stationary.taper.cov( x[1:3,],x[1:10,] , theta=1.5, Taper.args= list(k=2,theta=2.0,
                       dimension=2) )-> temp
# temp is now a tapered 3X10 cross covariance matrix in sparse format. 

 is.spam( temp)  # evaluates to TRUE

# should be identical to
# the direct matrix product

 temp2<- Exp.cov( x[1:3,],x[1:10,], theta=1.5) * Wendland(rdist(x[1:3,],x[1:10,]), 
                      theta= 2.0, k=2, dimension=2)
 test.for.zero(  as.matrix(temp), temp2)

# Testing that the "V" option works as advertized ...
x1<- x[1:20,]
x2<- x[1:10,]

V<- matrix( c(2,1,0,4), 2,2)
Vi<- solve( V)

u1<- t(Vi%*% t(x1))
u2<- t(Vi%*% t(x2))

look<- exp(-1*rdist(u1,u2))
look2<- stationary.cov( x1,x2, V= V)
test.for.zero( look, look2)


# Here is an example of how the cross covariance multiply works
# and lots of options on the arguments


 Ctest<- rnorm(10)
 
 temp<- stationary.cov( x,x[1:10,], C= Ctest, 
        Covariance= "Wendland", 
            k=2, dimension=2, theta=1.5 )

# do multiply explicitly

 temp2<- stationary.cov( x,x[1:10,],
        Covariance= "Wendland",
            k=2, dimension=2, theta=1.5 )%*% Ctest

 test.for.zero( temp, temp2)


# use the tapered stationary version 
# cov.args is part of the argument list passed to stationary.taper.cov
# within Krig. 
# This example needs the spam package.
# 

Krig(x,y, cov.function = "stationary.taper.cov", theta=1.5,
      cov.args= list(Taper.args= list(k=2, dimension=2,theta=2.0) )
           ) -> out2 
# NOTE: Wendland is the default taper here.

# BTW  this is very similar to 
Krig(x,y, theta= 1.5)-> out